function [correl_dX_dY,covar_dX_dY,std_dX,std_dY] = CalcTimeSeriesCorr2(X,Y,prob_X,Phi_X,n,Phi_X_n)
% computes corr(X(t+m)-X(t),Y(t+m)-Y(t))
    
    if nargin<=5
        Phi_X_n = mpower2(Phi_X,n);
    end
    
    mean_X = sum(prob_X(:).*X(:));
    var_X = sum(prob_X(:).*(X(:).^2)) - mean_X^2;
    mean_Y = sum(prob_X(:).*Y(:));
    var_Y = sum(prob_X(:).*(Y(:).^2)) - mean_Y^2;
    
    condE_X_n = Phi_X_n*X(:);
    var_dX = 2*var_X - 2*(sum(prob_X(:).*X(:).*condE_X_n) - mean_X^2);
    std_dX = sqrt(var_dX);
    
    condE_Y_n = Phi_X_n*Y(:);
    var_dY = 2*var_Y - 2*(sum(prob_X(:).*Y(:).*condE_Y_n) - mean_Y^2);
    std_dY = sqrt(var_dY);
    
    covar_dX_dY = sum(prob_X(:).*(2*X(:).*Y(:) - X(:).*condE_Y_n - Y(:).*condE_X_n(:)));
    
    correl_dX_dY = covar_dX_dY/(std_dX*std_dY);


end